

In [1]:

    
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt

from sympy import *
from mpmath import degrees, radians

init_printing()

This derivation is based on the USGS Professional Paper 1395: Map Projections — A Working Manual, ch 25, pp 195 – 199.



In [2]:

    
radius = Symbol('R')
semimajor, ellipticity = symbols('a e', positive=True)
central_lat, central_lon = symbols('phi_1 lambda_0')
point_lat, point_lon = symbols('phi lambda')

radius, semimajor, ellipticity, \
central_lat, central_lon, \
point_lat, point_lon









    Out[2]:





$$\left ( R, \quad a, \quad e, \quad \phi_{1}, \quad \lambda_{0}, \quad \phi, \quad \lambda\right )$$



In [3]:

    
def subdict(eq, vals):
    for k, v in vals.items():
        eq = eq.subs(k, v)
    return eq

Formulas for the Sphere



In [4]:

    
# Eq. 5-3
c = acos(sin(central_lat) * sin(point_lat) +
         cos(central_lat) * cos(point_lat) * cos(central_lon - point_lon))
c









    Out[4]:





$$\operatorname{acos}{\left (\sin{\left (\phi \right )} \sin{\left (\phi_{1} \right )} + \cos{\left (\phi \right )} \cos{\left (\phi_{1} \right )} \cos{\left (\lambda - \lambda_{0} \right )} \right )}$$



In [5]:

    
# Eq. 25-2
kp = c / sin(c)
kp









    Out[5]:





$$\frac{\operatorname{acos}{\left (\sin{\left (\phi \right )} \sin{\left (\phi_{1} \right )} + \cos{\left (\phi \right )} \cos{\left (\phi_{1} \right )} \cos{\left (\lambda - \lambda_{0} \right )} \right )}}{\sqrt{- \left(\sin{\left (\phi \right )} \sin{\left (\phi_{1} \right )} + \cos{\left (\phi \right )} \cos{\left (\phi_{1} \right )} \cos{\left (\lambda - \lambda_{0} \right )}\right)^{2} + 1}}$$



In [6]:

    
xs = radius * kp * cos(point_lat) * sin(point_lon - central_lon)
xs









    Out[6]:





$$\frac{R \sin{\left (\lambda - \lambda_{0} \right )} \cos{\left (\phi \right )} \operatorname{acos}{\left (\sin{\left (\phi \right )} \sin{\left (\phi_{1} \right )} + \cos{\left (\phi \right )} \cos{\left (\phi_{1} \right )} \cos{\left (\lambda - \lambda_{0} \right )} \right )}}{\sqrt{- \left(\sin{\left (\phi \right )} \sin{\left (\phi_{1} \right )} + \cos{\left (\phi \right )} \cos{\left (\phi_{1} \right )} \cos{\left (\lambda - \lambda_{0} \right )}\right)^{2} + 1}}$$



In [7]:

    
ys = radius * kp * (cos(central_lat) * sin(point_lat) -
                    sin(central_lat) * cos(point_lat) * cos(point_lon - central_lon))
ys









    Out[7]:





$$\frac{R \left(\sin{\left (\phi \right )} \cos{\left (\phi_{1} \right )} - \sin{\left (\phi_{1} \right )} \cos{\left (\phi \right )} \cos{\left (\lambda - \lambda_{0} \right )}\right) \operatorname{acos}{\left (\sin{\left (\phi \right )} \sin{\left (\phi_{1} \right )} + \cos{\left (\phi \right )} \cos{\left (\phi_{1} \right )} \cos{\left (\lambda - \lambda_{0} \right )} \right )}}{\sqrt{- \left(\sin{\left (\phi \right )} \sin{\left (\phi_{1} \right )} + \cos{\left (\phi \right )} \cos{\left (\phi_{1} \right )} \cos{\left (\lambda - \lambda_{0} \right )}\right)^{2} + 1}}$$

Formulas for the Ellipsoid

Unfortunately, the paper does not state full exact equations for the ellipsoid, but only approximations. I don't quite know how to write them in SymPy yet.

North Polar Aspect

The equation for $M$ is given as Eq. 3-21 on pg 197:

\begin{align} M = a [ & (1 - e^2 / 4 - 3e^4 / 64 - 5e^6 / 256 - \dots) \phi \\ - & (3e^2 / 8 + 3e^4 / 32 + 45e^6 / 1024 + \dots) \sin(2 \phi) \\ + & (15e^4 / 256 + 45e^6 / 1024 + \dots) \sin(4 \phi) \\ - & (35e^6 / 3072 + ...) \sin(6 \phi) \\ + & \dots ] \end{align}

The numerator in the first term can be written easily as: $$ \sum_{n=1, 2, \dots} (2n-1) e^{2n} $$ but the denominator is: $$ 4, 64, 256, \dots = 2^2, 2^6, 2^8, \dots $$ which doesn't seem to follow a sane pattern, and the second, third, fourth, etc. terms are worse. Possibly, some trig rules are reducing things, but the pattern is not particularly evident the way they've written it. I also haven't yet read chapter 3, so perhaps this is all explained there.

Let's try out the reduced form anyway, assuming they kept only the relevant terms...



In [8]:

    
M = semimajor * (
    (1 - ellipticity**2 / 4 - 3 * ellipticity**4 / 64 - 5 * ellipticity**6 / 256) * point_lat -
    (3 * ellipticity**2 / 8 + 3 * ellipticity**4 / 32 + 45 * ellipticity**6 / 1024) * sin(2 * point_lat) +
    (15 * ellipticity**4 / 256 + 45 * ellipticity**6 / 1024) * sin(4 * point_lat) -
    (35 * ellipticity**6 / 3072) * sin(6 * point_lat)
)
M









    Out[8]:





$$a \left(- \frac{35 e^{6}}{3072} \sin{\left (6 \phi \right )} + \phi \left(- \frac{5 e^{6}}{256} - \frac{3 e^{4}}{64} - \frac{e^{2}}{4} + 1\right) + \left(\frac{45 e^{6}}{1024} + \frac{15 e^{4}}{256}\right) \sin{\left (4 \phi \right )} - \left(\frac{45 e^{6}}{1024} + \frac{3 e^{4}}{32} + \frac{3 e^{2}}{8}\right) \sin{\left (2 \phi \right )}\right)$$



In [9]:

    
Mp_np = M.subs(point_lat, radians(90))
Mp_np









    Out[9]:





$$a \left(- 0.0306796157577128 e^{6} - 0.0736310778185108 e^{4} - 0.392699081698724 e^{2} + 1.5707963267949\right)$$



In [10]:

    
rho_np = Mp_np - M
rho_np









    Out[10]:





$$a \left(- 0.0306796157577128 e^{6} - 0.0736310778185108 e^{4} - 0.392699081698724 e^{2} + 1.5707963267949\right) - a \left(- \frac{35 e^{6}}{3072} \sin{\left (6 \phi \right )} + \phi \left(- \frac{5 e^{6}}{256} - \frac{3 e^{4}}{64} - \frac{e^{2}}{4} + 1\right) + \left(\frac{45 e^{6}}{1024} + \frac{15 e^{4}}{256}\right) \sin{\left (4 \phi \right )} - \left(\frac{45 e^{6}}{1024} + \frac{3 e^{4}}{32} + \frac{3 e^{2}}{8}\right) \sin{\left (2 \phi \right )}\right)$$



In [11]:

    
xnp = rho_np * sin(point_lon - central_lon)
ynp = -rho_np * cos(point_lon - central_lon)

South Polar Aspect



In [12]:

    
Mp_sp = M.subs(point_lat, radians(-90))
Mp_sp









    Out[12]:





$$a \left(0.0306796157577128 e^{6} + 0.0736310778185108 e^{4} + 0.392699081698724 e^{2} - 1.5707963267949\right)$$



In [13]:

    
rho_sp = Mp_sp + M
rho_sp









    Out[13]:





$$a \left(0.0306796157577128 e^{6} + 0.0736310778185108 e^{4} + 0.392699081698724 e^{2} - 1.5707963267949\right) + a \left(- \frac{35 e^{6}}{3072} \sin{\left (6 \phi \right )} + \phi \left(- \frac{5 e^{6}}{256} - \frac{3 e^{4}}{64} - \frac{e^{2}}{4} + 1\right) + \left(\frac{45 e^{6}}{1024} + \frac{15 e^{4}}{256}\right) \sin{\left (4 \phi \right )} - \left(\frac{45 e^{6}}{1024} + \frac{3 e^{4}}{32} + \frac{3 e^{2}}{8}\right) \sin{\left (2 \phi \right )}\right)$$



In [14]:

    
xsp = xnp
ysp = rho_sp * cos(point_lon - central_lon)

Guam Projection



In [15]:

    
M_guam = M
M1_guam = M_guam.subs(point_lat, central_lat)
M1_guam









    Out[15]:





$$a \left(- \frac{35 e^{6}}{3072} \sin{\left (6 \phi_{1} \right )} + \phi_{1} \left(- \frac{5 e^{6}}{256} - \frac{3 e^{4}}{64} - \frac{e^{2}}{4} + 1\right) + \left(\frac{45 e^{6}}{1024} + \frac{15 e^{4}}{256}\right) \sin{\left (4 \phi_{1} \right )} - \left(\frac{45 e^{6}}{1024} + \frac{3 e^{4}}{32} + \frac{3 e^{2}}{8}\right) \sin{\left (2 \phi_{1} \right )}\right)$$



In [16]:

    
x = (semimajor * (point_lon - central_lon) * cos(point_lat) /
     (1 - ellipticity**2 * sin(point_lat)**2)**0.5)
y = (M_guam - M1_guam + x**2 * tan(point_lat) *
     (1 - ellipticity**2 * sin(point_lat)**2)**0.5 / (2 * semimajor))

Cartopy Boundaries



In [17]:

    
# Default Globe in Cartopy
WGS84_semimajor = Rational(6378137)
WGS84_flattening = Rational(1, 298.257223563)
# Plotting fails if I leave this bit as exact, thinking there are complex values,
# even though the input to sqrt is positive, and it's supposed to return the positive root.
WGS84_ellipticity = sqrt(2 * WGS84_flattening - WGS84_flattening**2).evalf()

# Default parameters in Cartopy
lat0_default = radians(0)
lon0_default = radians(0)

Sphere



In [18]:

    
# From numerical examples, pg 337
radius_example = 3.0
lat0_example = 40
lon0_example = -100

sphere_vars = {radius: radius_example,
               central_lat: radians(lat0_example),
               central_lon: radians(lon0_example)}

I chose the following boundary locations since they work for other projections, but I'm not certain they'll work here.



In [19]:

    
lon_min_sphere = radians(lon0_example - 180)
lon_max_sphere = radians(lon0_example + 180)
lat_min_sphere = radians(lat0_example - 90)
lat_max_sphere = radians(lat0_example + 90)

example_xs = subdict(xs, sphere_vars)
eq1x = example_xs.subs(point_lat, lat_max_sphere)
eq2x = example_xs.subs(point_lon, lon_max_sphere)
eq3x = example_xs.subs(point_lat, lat_min_sphere)
eq4x = example_xs.subs(point_lon, lon_min_sphere)

example_ys = subdict(ys, sphere_vars)
eq1y = example_ys.subs(point_lat, lat_max_sphere)
eq2y = example_ys.subs(point_lon, lon_max_sphere)
eq3y = example_ys.subs(point_lat, lat_min_sphere)
eq4y = example_ys.subs(point_lon, lon_min_sphere)

plotting.plot_parametric((eq1x, eq1y, (point_lon, lon_min_sphere, lon_max_sphere)),
                         (eq2x, eq2y, (point_lat, lat_max_sphere, lat_min_sphere)),
                         (eq3x, eq3y, (point_lon, lon_max_sphere, lon_min_sphere)),
                         (eq4x, eq4y, (point_lat, lat_min_sphere, lat_max_sphere)),
                         title='Guessed Sphere Boundaries')









    












    Out[19]:





<sympy.plotting.plot.Plot at 0x7f0248b287f0>



In [20]:

    
lat_magnitude = subdict((point_lat - central_lat) / radians(90), sphere_vars)
lon_magnitude = subdict((point_lon - central_lon) / radians(180), sphere_vars)

plotting.plot3d_parametric_surface(example_xs, example_ys, lat_magnitude,
                                   (point_lat, lat_min_sphere, lat_max_sphere),
                                   (point_lon, lon_min_sphere, lon_max_sphere),
                                   xlabel='X', ylabel='Y', title='Sphere - Latitude Dependence')

plotting.plot3d_parametric_surface(example_xs, example_ys, lon_magnitude,
                                   (point_lat, lat_min_sphere, lat_max_sphere),
                                   (point_lon, lon_min_sphere, lon_max_sphere),
                                   xlabel='X', ylabel='Y', title='Sphere - Longitude Dependence')









    












    












    Out[20]:





<sympy.plotting.plot.Plot at 0x7f0248aec588>

The longitude and latitude dependence is pretty complicated (worse than the above seems to show). Maybe we can find the boundaries using the derivative, but I'm not sure.



In [21]:

    
xs_deriv_lat = diff(xs, point_lat)
xs_deriv_lon = diff(xs, point_lon)
ys_deriv_lat = diff(ys, point_lat)
ys_deriv_lon = diff(ys, point_lon)



In [22]:

    
xsd_lat = subdict(xs_deriv_lat, sphere_vars)
xsd_lon = subdict(xs_deriv_lon, sphere_vars)
ysd_lat = subdict(ys_deriv_lat, sphere_vars)
ysd_lon = subdict(ys_deriv_lon, sphere_vars)



In [23]:

    
# Solving them together uses 14+ GiB and still doesn't complete!
# solve((xsd_lat, xsd_lon), (point_lat, point_lon))

# extreme_lat = solve(xsd_lat.subs(point_lon, 0), point_lat)

The above solve does not complete in either manner. So I'm going to bench this attempt for the time being.

Ellipsoid



In [24]:

    
ellipsoid_vars = {semimajor: WGS84_semimajor,
                  ellipticity: WGS84_ellipticity,
                  central_lat: radians(lat0_default),
                  central_lon: radians(lon0_default)}

North Polar Aspect



In [25]:

    
lat0_np = 90
lon0_np = 0

np_vars = {semimajor: WGS84_semimajor,
           ellipticity: WGS84_ellipticity,
           central_lat: radians(lat0_np),
           central_lon: radians(lon0_np)}

lon_min_np = radians(lon0_np - 180)
lon_max_np = radians(lon0_np + 180)
lat_min_np = radians(lat0_np - 180)
lat_max_np = radians(lat0_np)

default_xnp = subdict(xnp, np_vars)
default_ynp = subdict(ynp, np_vars)

lat_magnitude = subdict(point_lat / radians(90), np_vars)
lon_magnitude = subdict(point_lon / radians(180), np_vars)

plotting.plot3d_parametric_surface(default_xnp, default_ynp, lat_magnitude,
                                   (point_lat, lat_min_np, lat_max_np),
                                   (point_lon, lon_min_np, lon_max_np),
                                   xlabel='X', ylabel='Y',
                                   title='Ellipsoid - North Polar Aspect - Latitude Dependence')

plotting.plot3d_parametric_surface(default_xnp, default_ynp, lon_magnitude,
                                   (point_lat, lat_min_np, lat_max_np),
                                   (point_lon, lon_min_np, lon_max_np),
                                   xlabel='X', ylabel='Y',
                                   title='Ellipsoid - North Polar Aspect - Longitude Dependence')









    












    












    Out[25]:





<sympy.plotting.plot.Plot at 0x7f024dd25f60>

From the above, it appears the boundary occurs on the minimum latitude, at least for the North Polar Aspect.

South Polar Aspect



In [26]:

    
lat0_sp = -90
lon0_sp = 0

sp_vars = {semimajor: WGS84_semimajor,
           ellipticity: WGS84_ellipticity,
           central_lat: radians(lat0_sp),
           central_lon: radians(lon0_sp)}

lon_min_sp = radians(lon0_sp - 180)
lon_max_sp = radians(lon0_sp + 180)
lat_min_sp = radians(lat0_sp)
lat_max_sp = radians(lat0_sp + 180)

default_xsp = subdict(xsp, sp_vars)
default_ysp = subdict(ysp, sp_vars)

lat_magnitude = subdict(point_lat / radians(90), sp_vars)
lon_magnitude = subdict(point_lon / radians(180), sp_vars)

plotting.plot3d_parametric_surface(default_xsp, default_ysp, lat_magnitude,
                                   (point_lat, lat_min_sp, lat_max_sp),
                                   (point_lon, lon_min_sp, lon_max_sp),
                                   xlabel='X', ylabel='Y',
                                   title='Ellipsoid - South Polar Aspect - Latitude Dependence')

plotting.plot3d_parametric_surface(default_xsp, default_ysp, lon_magnitude,
                                   (point_lat, lat_min_sp, lat_max_sp),
                                   (point_lon, lon_min_sp, lon_max_sp),
                                   xlabel='X', ylabel='Y',
                                   title='Ellipsoid - South Polar Aspect - Longitude Dependence')









    












    












    Out[26]:





<sympy.plotting.plot.Plot at 0x7f02486b28d0>

Again, the boundary occurs at the minimum latitude.

Oblique Aspect

After some thought, it seems that the boundary should be at the greatest distance away from the central point. This makes sense given the properties of this projection. The projection plots azimuth after all. This even explains the "strangeness" for the spherical case; because the central point is not at the poles, the radial distance must be greater than zero. Just blindly taking negative and positive latitudes and longitudes will not cover the whole azimuth, and sometimes produces odd wraparound.

But this actually is good news! It means the boundary is simply half a circumference. Obviously, the meaning is slightly different for an ellipsoid, but it's certainly a lot easier to calculate.